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O Abstract 
<U 

^ Elementary smooth functions (beyond contact) are employed to construct pair correlation func- 

tions that mimic jammed disordered sphere packings. Using the 52-invariant optimization method 

^ of Torquato and Stillinger [J. Phys. Chem. B 106, 8354, 2002], parameters in these functions are 

g optimized under necessary realizability conditions to maximize the packing fraction (f) and average 

number of contacts per sphere Z. A pair correlation function that incorporates the salient features 

O of a disordered packing and that is smooth beyond contact is shown to permit a ^ of 0.6850: 

this value represents a 45% reduction in the difference between the maximum for congruent hard 

^ spheres in three dimensions, 7r/vl8 0.7405, and 0.64, the approximate fraction associated with 

O . 

maximally random jammed (MRJ) packings in three dimensions. We show that, surprisingly, the 

m 

CO continued addition of elementary functions consisting of smooth sinusoids decaying as r ^ permits 

^ packing fractions approaching 7r/vl8. A translational order metric is used to discriminate between 
00 

O degrees of order in the packings presented. We find that to achieve higher packing fractions, the 

. !^ degree of order must increase, which is consistent with the results of a previous study [Torquato 

?H et al., Phys. Rev. Lett. 84, 2064, 2000]. 

PACS numbers: 
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I. INTRODUCTION AND BACKGROUND 



Packing problems address the various arrangements of a set (finite or infinite) of non- 
overlapping objects in a space of given dimension [U 121 E]- Often, one seeks to arrange the 
objects in such a way as to optimize a certain statistical or bulk property, e.g. the number 
density p of objects (or equivalently the packing fraction 0, the fraction of space covered by 
the objects' interiors). This paper is concerned with optimizing the packing fraction cj) and 
average number of contacts per sphere (average kissing number) Z for packings of congruent 
spheres in three-dimensional Euclidean space while maintaining an assumed functional form 
for the pair correlation function (defined below) of a "random" packing. 

Random packings of three-dimensional hard spheres have been studied by scientists to 
better understand everything from heterogeneous materials to liquids to granular media (like 
sand) to living cells [Sj HJ |5l [6l [?j . In parallel to the concept of a maximum packing fraction 
for periodic (crystalline) packings, it had been assumed that a (different) maximum packing 
fraction could be defined for random packings, referred to as "random close packing" (RCP) 
[HI E] • However, while there is a proved maximum packing fraction for hard-sphere periodic 
packings achieved via an FCC lattice or one of its stacking variants [2j , the RCP state has 
been shown to be ill-defined [lOl [TT] . 

Experimental packings of oiled steel ball bearings originally led to the idea that me- 
chanically stable random packings of identical spheres could not exhibit packing fractions 
exceeding 0.64 or declining past 0.60 [HI [12] • Mathematically constructed models P3] and 
early computer simulations [H] seemed to support these conclusions, though later work 
demonstrated that the limiting packing fractions obtained were highly dependent on the 
packing methods [3l[T0]. These methods included, for example, lightly vibrating a container 
filled with spheres in either horizontal or vertical motions, rolling spheres one by one into a 
container [T^j, and simulating the compression of a hard-sphere gas [TU| [T7]. 

Changing the method of packing via molecular dynamics simulations showed that den- 
sities past 0.64 are realizable [ini [iHl HE]- However, it has yet to be demonstrated from 
a theoretical basis, i.e., without resort to experiment or computer simulation, that indeed 
there is no maximum density limit for "random" packings apart from the proved limit for 
periodic packings. This paper provides such a theoretical basis by extending previous op- 
timization studies of (yf2-invariant processes (defined below) pUl [21] to a broader class of 
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disordered packings. 

Following previous work by two of us, a statistically homogeneous and isotropic packing 
is defined to be disordered if its pair correlation function (72 (r) in (i-dimensional Euclidean 
space M'^ decays to unity faster than r~'^~'^ for some e > [21j. Each packing corresponds 
to a unique g2{r), a function proportional to the probability density of finding a separation 
r between any two sphere centers and normalized such that it takes the value of unity 
when no spatial correlations are present. The precise definition of "disordered" via the pair 
correlation function takes the place of the imprecise term "random" hereafter. 

The essential ideas behind our approach were actually laid out in our earlier work [201 EI] ■ 
In Ref. |20], the main objective was to study disordered packings in which short-range order 
was controlled using so-called (yf2-invariant processes. A g2-invariant process is one in which 
a given pair correlation function g2{T) remains invariant for all r as packing fraction varies 
over the range of densities 

< </) < (1) 

The terminal packing fraction 0=,, is the maximum achievable for the (yf2-invariant process 
subject to satisfaction of the nonnegativity of g2{f) and the structure factor S{k), i.e., 

g2{r) > V r > 0, (2) 

S{k) > V > 0. (3) 

For a statistically homogeneous and isotropic packing at number density p, S{k) is related 
to the Fourier transform of the total correlation function h{r) = g2{r) — 1, by 

Sik) = l + ph{k), (4) 

where h{k) represents the Fourier transform of h{r). In three dimensions, this can be written 

as 

S{k) = 1 + 47rp r '^^^h{r)dr. (5) 
Jo k 

The optimization procedure described above was formulated for a hard-sphere packing, 
which requires the additional condition on (72 that respects the nonoverlap constraint, i.e., 

g^(r) = for 0<r < D. (6) 

When there exist sphere packings with a (72 satisfying conditions ([2]), (|3]), and ^ for 
in the interval [0,0^,], then a lower bound on the maximal packing fraction is given by 



''max 



> 0*. It is noteworthy that this optimization problem for sphere packings is an 
infinite-dimensional linear program, which is the dual of the primal linear program devised 
by Cohn and Elkies [22] to obtain upper bounds on the maximal packing fraction [2T]. We 
will comment further on this connection in the conclusions. Finally, we note here that the 
results of the optimization of the pair correlation function given in Ref. [20] has found 
application in describing small-scale convective structural features of the solar surface [23] . 

The nonnegativity conditions ^ and (|3| are necessary, but generally not sufficient, for 
a pair correlation function at a given density to be realizable by a point process [21]. A 
third condition, obtained by Yamada [25] and not included in the optimization procedure 
described above, constrains cr'^{A) = {{N{A) — {N{A)))'^), which is the variance in the 
number of points N{A) contained within a window A G M''; 



'iA)>e{A)[i-e{A)], (7) 



where 0{A) is the fractional part of the expected number of points contained within the win- 
dow. The number variance associated with a spherical window of radius R for a statistically 
homogeneous point process in d dimensions can be written as follows [26] : 

a^R) = pvi{R) [1 + p [ /i(r)4"*(r; R)dr] > 9{R) [l - d{R)\ , (8) 

where vi{R) is the volume of the window and a'^^^{r;R) the intersection volume of two 
windows of radius R (whose centers are separated by r) divided by vi{R). This additional 
condition ([S]) on the pair correlation function is always satisfied for statistically homogeneous 
and isotropic packings with sufficiently large windows in dimensions greater than 1 [21j. In 
all cases that have been studied, this condition is satisfied for all R if the first two conditions 
are satisfied [211127]. 

While conditions g, Q, and Q are necessary for the realizability of point processes, 
along with incorporation in (72 (r) of the core exclusion feature, they appear to be rather 
strong conditions for realizability of sphere packings, especially as the space dimension in- 
creases [21]. For example, a method to construct disordered packing configurations that 
realize test 5'2's meeting the conditions and incorporating the features of core exclusion and 
contact pairs [Eqs. Q and (|lT|)] has been successful [2HI EH] in two and three dimensions. 



No example in three dimensions or greater of an unrealizable g2 satisfying the conditions 
and incorporating the core exclusion feature is currently known. 



In Ref. [20], a five-parameter test family of g2^s incorporating features of core exclusion, 
contact pairs, and damped oscillatory short-range order beyond contact [Eqs. ([9]), (11), 
and (12)] had been considered. The problem of finding the terminal packing fraction 0^, was 
posed as an optimization problem: maximize over the set of parameters subject to the first 
two realizability conditions (the third condition due to Yamada was not relevant). In this 
work, we consider a broader family of smooth g2 test functions corresponding to disordered 
packings [30] and satisfying all three aforementioned conditions. 

To demonstrate the absence of a theoretical upper limit on disordered packings, we show 
that terms decaying as r~^, representative of a feature prominent in the pair correlation 
functions of maximally random jammed (MRJ) packings [31], allow for increased packing 
fraction for pair correlation functions satisfying the three conditions and incorporating the 
aforementioned features. A simple 11-parameter form consisting of the initial five-parameter 
form plus two sinusoids decaying as permits a packing fraction of 0.6850. Using a trans- 
lational order metric, we show that the pair correlation function with the highest packing 
fraction also exhibits the highest degree of order, which is consistent with the conclusions of 
a previous work [lOj. Additionally we show the surprising result that the continued addition 
of terms decaying as allows for packing fractions up to vr/v^TS, indicating that, if the 
packings are realizable, the progression of disordered packings up to the maximum is a 
continuum, dependent only on the form and parameters of the functions employed. A qual- 
itative description of a realizable disordered packing with smooth (72 (r) and approaching 
tt/vTS is provided in Section IV A. 



II. OPTIMIZATION OF 52-INVARIANT PROCESSES 

We begin by revisiting the optimization problem first examined in Ref. [20] • We employ 
a more comprehensive search using simulated annealing to optimize the five parameters 
of the family of (72 's presented in Ref. [20] and find a higher terminal packing fraction 
0* = 0.64268. The three functions that comprise the five parameter family, gi{r), gn{r) 
and giiiir), capture the most salient properties of a disordered packing, including that the 
average number of spheres in contact is Z and that no sphere centers may approach closer 
than a distance of one sphere diameter. 
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A Heaviside step function represents the spheres' hard core exclusion, 

y/ = e(r-l), (9) 

where we set the diameter of the spheres to be unity. The Heaviside step function Q{x) is 
defined piecewise as 

{0, X <0 
(10) 
1, x>0. 

A Dirac delta function represents pair contacts, 

9n = ^^Sir-l), (11) 

with Z the average number of contacts per sphere (average kissing number). An exponen- 
tially decaying sinusoid provides short-range oscillatory motion about unity: 

gni = ^ exp (-Sir) sin(Cir + D,)e{r - 1), (12) 

with parameters Ai, Bi, Ci, Di. The total pair correlation function g2{T) is then 

g2{r) = gi{r) + gn{r) + giii{r). (13) 

Constraining i?i > 0, Z > ensures a physical configuration, while constraining Ci > 0, 
< < TT eliminates function duplicates without additionally constraining the range of 
the functional form. 



A. Maximizing packing fraction 

Packing fraction is maximized using the 5f2-iiivariant method in 20,000 independent runs 

of over 10,000 iterations each. Initial parameters arc confined to the following bounds, and 
selected randomly before each run with exponentially decreasing probability from zero: 

-50 <Ai< 100 
< Si < 10 
< Ci < 50 

< L>i < TT 

< Z < 13 
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Parameters are allowed to range outside of bounds, but in no cases of the 20,000 did this 
occur. 

Using Gj{k), Gjjik), Guiik) to represent 1/p times the second term on the right hand 
side of relation ([s]) with h{r) = gi{r) + gn{r) + giii{r) — 1, the structure factor for the 
functions becomes 

S^k) = 1 + + p{Gi{k) + Gui{k)), (14) 

with Gjj{k) = Zsm{k)/pk. The exact analytical forms for Gj, Gjj, and Gju are included 
in Appendix B. 

The method to maximize relies upon a simple principal. If S{K) = at some point 
K and Gi{K) + Giii{K) < 0, while S{k) > for all other points k, then is at a global 
maximum for the given five parameters, i.e., = 0*. Hence to maximize 0*, S{k) is 



analytically calculated from the pair correlation function in accordance with Eq. (14), and 
for each random step along one of the five parameters, if possible is chosen such that the 
structure factor is in accordance with this principal. If obeying the structure factor condition 
is not mathematically possible for the parameter set, if the pair correlation function g2{r) 
is not greater than or equal to zero for all r, or if the maximum for the set does not meet 
the standard temperature-dependent simulated annealing condition for accepting a move, 
the random step is rejected and a new step chosen. 

Figures [T] and [2] present pair correlation functions and corresponding structure factors that 
yield terminal packing fractions 0=,, = 0.64268 and 0=,, = 0.64050, respectively. Parameters 
for high packing-fraction results were similar in period Ci, phase Di, and average kissing 
number Z, but varied in amplitude Ai and damping factor Bi. For example, the g2{f') 
with the highest two terminal packing fractions, 0* = 0.64268 and 0* = 0.64050, exhibited 
Ai = 3.0996 and Ai = 6.1707, with Bi = 0.58091 and Bi = 1.2090, respectively. Each 
g2{r) exhibits a minimum near r = 1.35, with the minimum for the highest packing fraction 
equal to zero, and a maximum at r = 1 of about 2.7, suggesting that these traits, along 
with a period 2tt/Gi of about 1.2 and a phase Di of about 0.50 are important in obtaining 
the maximum packing fraction for this functional form. Further analysis indicates that the 
period and existence of a deep minimum within r 1.5 remain important in maximizing 
packing fraction when other elements are added to this functional form. 
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FIG. 1: Top: The pair correlation function g2{r) for (j)^ = 0.64268, Ai = 3.0996, Bi = 0.58091, 
Ci = 7.54069, Di = 0.45970, Z = 5.0633. Bottom: The corresponding structure factor ^(A;). 

B. Maximizing kissing number 



Maximum packing-fraction g2^s for this functional form do not correspond to g2^s that 
maximize average kissing number, as is the general case for sphere packings in many dimen- 
sions |27]. Moreover, though the average kissing number for the highest possible packing 
fraction (FCC lattice) is 12 [32j, the proved maximum possible, Z may only obtain the value 
of 9.5401 for this form. As will be seen later, as additional elements and parameters are 
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FIG. 2: Top: The pair correlation function g2{r) for (j)^ = 0.64050, Ai = 6.1707, Bi = 1.2090, 
Ci = 7.6011, Di = 0.54981, Z = 5.1593. Bottom: The corresponding structure factor S{k). 



added to this form and optimized for maximum packing fraction, average kissing number 
increases substantially. It is of interest therefore to maximize Z with packing fraction as 
a parameter. The method employed is the same as before: simulated annealing in 20,000 
independent runs. 

Figure [3] shows the pair correlation function for the optimal parameters Ai, Bi, Ci, Di, 
and 0, with optimized Z = 9.5401. Most notably, the packing fraction associated with 
the maximum average kissing number configuration is 0.631, nearly equal to the maximum 
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FIG. 3: The pair correlation function g2{r) for Z = 9.532, Ai = -29.02, Si = 2.735, Ci = 6.0537, 
Di = 0.47, (/) = 0.631 

packing fraction achieved (0=,, = 0.64268) for this functional form. The fact that many more 
spheres are in contact (relative to the case in which the packing fraction is maximized, where 
Z = 5.0633), could reflect the presence of large clusters in the packing or even a sample- 
spanning cluster [23] • Also of note is that the first minimum, while still equal to zero, has 
moved inwards toward r = 1, and that the value of the pair correlation function at r = 1+, 
g2{^~^) = 0.5514, is much smaller than before, where the notation 1"*" refers to the right hand 
limit of g2{r) as r ^ 1. 

These features must be present with an average kissing number near 12, as the integral 
of g2{r) from r = 1^ to a given R is related to the number of sphere centers most likely to 
be found in that range of distances, and for a realizable packing, this number cannot exceed 
12 for R near unity. Specifically, Z{R), the total expected number of sphere centers to be 
found within a (larger) sphere of radius R centered on an arbitrary sphere center within the 
packing, can be written 



where the average of Z{R) as R —>■ 1 from the right (1+) is equivalent to the average kissing 
number Z. 

Generally, Z{R) cannot exceed some ZmaxiR), where Zmax{R) represents the maximum 




(15) 
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number of sphere centers that can be placed within a (larger) sphere of radius R centered 
on an arbitrary sphere center. It is clear from geometric considerations that Z^axiR) is 
a piecewise continuous nowhere decreasing function of R. As the maximum number of 
congruent spheres that can be placed around a congruent central sphere is 12, Zmax{R) is 
12 on the interval {R : (1,1 + a)} for some small deterministic parameter a > 0, where a 
is the distance from unity to the first discontinuity in Zmax{R)- Currently, a is not known 
rigorously, though its value is suspected to be about 0.045 [31]. That Z{R) must be less than 
or equal to Zmax{R) for all R is another necessary realizability condition for (72 's representing 
sphere packings. More will be said about this condition in the conclusions. 

It is also of note that the damping parameter Bi for the maximum average kissing number 
case is much larger than for the maximum packing fraction case, implying that in the max- 
imum average kissing number case, spatial correlation decreases substantially more quickly. 
The configuration represented by the more heavily damped pair correlation function can be 
interpreted as many groups of tightly packed spheres pushed together in a random fashion, 
thus exhibiting a high average kissing number but little to no spatial correlation at distances 
greater than several sphere diameters. 

III. INCREASING PACKING FRACTION 

Here we show that the addition past contact of sinusoids decaying (to unity) as T , a 
sinusoidally modulated form of the r~'^ decay present in pair correlation functions calculated 
from large-scale simulations (10® spheres) of MRJ-like packings |3T], to the initial five- 
parameter family of (72's allows for significant increase in terminal packing fraction above 
0.64268. In three dimensions, the inverse 4*'^ power is the smallest integer power to which a 
pair correlation function may decay while satisfying the r~^~^ condition for representing a 
disordered packing. However, aside from the clues that this mathematical attribute might 
provide, the structural origins of the decay currently remain a conceptual mystery. 

While the form of r~'^ decay allows increased packing fraction, other forms, including 
those representative of certain other features present in MRJ states, do not necessarily. 
In fact many other additions have been considered for this study, though from these no 
substantial increases in 0=,, above the 0.64268 achieved for the five-parameter form have been 
obtained. This is not to say that only sinusoids decaying as will increase the maximum 
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possible attainable packing fraction, just that the selection of elements is non-trivial. In 
particular, one feature present in MRJ states, a fractional power-law divergence near r = 1, 
resulted only in a reduced value for (p^. The addition of this feature is discussed in Appendix 
B. 

The method used to optimize parameters with additional forms included involves three 



steps. Structure factors are calculated analytically from Eq. (14) as before, but to accommo- 
date the increased processing time required to find the minima of S{k) due to the complexity 
of its analytic form, the number of points k at which S{k) is calculated is initially reduced. 
The first step then is to find several "rough" maximum density configurations, just as with 
the initial five-parameter configurations, using a reduced number of calculations with initial 
parameters selected as before in 10,000 independent runs. The second step is to improve 
upon these approximate maxima by increasing the number density of points k calculated in 
S{k) about its minimum, in 1,000 additional runs for each rough maximum 0. The results 
of the second step still do not yield an exact maximum and hence the third step is to fine- 
tune the maximum from the second step manually, ensuring that S{k) and g2{r) are indeed 
greater than zero for all r and k. 

Sinusoids decaying as are highly successful in increasing maximum packing fraction 
beyond 0.64268. Specifically, cosine functions of the form 

gjy{r) = ^cos{Br + C)e{r-l) (16) 

are employed. Adding two of these elements to the initial five-parameter form allows iden- 
tification of a function obeying all conditions with a maximum packing fraction of 0.6850, 
45% closer to the packing fraction of the FCC configuration than 0.640. Function (17) is 
the 11-parameter pair correlation functional form specified. 

^2(r) = -^Sir -1)+ + ^ exp {-Bir) sin(Cir + Di)^ e(r - 1) 

+ cos(52r + C2) + ^ cos{Bsr + Q)^ Q{r - 1), (17) 

or 

92ir) = gi{r) + gji{r) + gmir) + giva{,r) + givbir), (18) 
where the subscripts a and b in the last two terms on the right refer to the fourth and fifth 



terms in (17). The analytical structure factors for the pair correlation functions represented 
by function (17) calculated using relation ([s]) are given in Appendix A. 
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FIG. 4: The pair correlation function g2{r) for = 0.6850, Ai = 3.555, Bi = 0.7189, Ci = 7.5589, 
Di = 0.6247, Z = 10.4281, A2 = 0.2205, ^2 = 2.370, C2 = 0.000, A3 = 2.2822, S3 = 8.7373, 
C3 = 0.0423 

Figures ^ and (|5| show graphs of the pair correlation function and structure factor with 



the 11 parameters of function (17) optimized for a maximum packing fraction of 0.6850. 



The first minimum and height of g2{r) just past r = 1 in Fig. |4]are similar to those observed 
in Fig. |3j the plot of the pair correlation function for a maximized average kissing number 
using the initial five-parameter functional form. The period 27r/Ci, phase Di, and decay 
factor Bi of the exponentially decaying sinusoidal function are similar to those observed in 
the top plots of Figs. [T] and |2| the plots of the pair correlation functions for maximized 
packing fraction using the initial five-parameter functional form. The minimum and height 
just past r = 1 must be similar for physical meaningfulness due to high average kissing 
number Z = 10.4281, while the similarities in period, phase, and decay factor indicate that 
a higher degree of order is maintained farther from r = 1. 

The decay factor, Bi = 0.7189, for the 0.6850 maximum packing fraction function is in 
fact larger than the decay factor, Bi = 0.58091, for the 0.64268 packing fraction function, 
but fast decay in the former g2 function is avoided through the addition of the sinusoids. 

The average kissing number Z = 10.4281 for the 0.6850 maximum packing fraction func- 
tion is substantially higher than that obtained for the 0.6427 maximum packing fraction 
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FIG. 5: The structure factor S{k) for cj)^ = 0.6850, Ai = 3.555, Bi = 0.7189, Ci = 7.5589, 
Di = 0.6247, Z = 10.4281, A2 = 0.2205, B2 = 2.370, C2 = 0.000, A3 = 2.2822, S3 = 8.7373, 
C3 = 0.0423 

function. The addition of the elements allows not only for tighter packing of spheres 
(and hence a higher Z than previously possible), but also for slower decay than present 
with an exponentially decaying function alone. This implies (but does not prove) that more 
correlation at greater distances is necessary for higher packing fraction disordered packings: 
a hypothesis that will be supported further later in this paper by order metric calculations. 

It is of note that reducing the packing fraction without proportionally reducing the 
kissing number Z quickly violates the structure factor condition. This implies that kissing 
numbers as high as 10.42 cannot be maintained without proportionally high densities, which 
is in agreement with physical intuition: average number of spheres in contact cannot increase 
past a certain point without high enough packing fraction, though it is important to note 
that the converse of this statement is not true (if a configuration is not required to be 
jammed). 

Further supporting this notion is Fig. [6| a graph of the pair correlation function with 
10 parameters optimized to maximize with Z set to 6.672. With the four parameters of 



function (12) nearly the same as for the 0.6850 terminal packing fraction function but Z 



reduced substantially, only a lower terminal packing fraction is possible. The amplitudes of 
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FIG. 6: The pair correlation function 52 (r) for cj)^ = 0.6700, Ai = 3.3048, Bi = 0.71901, Ci = 
7.5115, Di = 0.6610, Z = 6.672, A2 = 0.0795, B2 = 2.390, C2 = 0.000, A3 = 0.312, B3 = 8.389, 
C3 = 0.000 



the cosine elements are also substantially reduced in Fig. |6| implying that the increase in 
packing fraction from 0.6700 to 0.6850 necessitates both correlation at greater distances and 
a substantially higher average kissing number. As will be quantified later in this paper, it 
follows that more order overall is required to attain the packing fraction 0.6850; nonetheless, 
the smooth correlation function past contact and swift decay to unity indicates that the 
configuration remains disordered as defined. 

As before with the five-parameter functional form, it is of interest to maximize kissing 



number Z for the 11-parameter functional form (17). We find that the average kissing 



number cannot increase much above 10.4281. For this functional form and with realizability 
conditions satisfied, we see that if one of either the packing fraction or kissing number is 
near its maximum value, then the other must also be near its maximum. 
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IV. FROM DISORDERED TO CRYSTALLINE 



A. Order metrics 

A disordered packing as defined clearly lacks the degree of order present in a periodic 
packing, but the degree of order within a disordered system still varies, and may be quan- 
tified. For cases of hard-sphere packings, many different methods of quantifying the order 
in a configuration are possible. The majority of these methods advocate the use of scalar 
statistical measures. Some of these assume that the most ordered system is the appropriate 
maximally dense packing [31 HDl [ISl [121 ES]' while others do not presuppose a reference 
crystal state P [TBI CHI ESI |36] . 

For our purposes, the translational order metric introduced by Truskett et al. [18] is 
convenient because it is given in terms of the total correlation function: 



1 



^/ \h{xp'"')\dx. (19) 



Xc - P^''^D Jpi/3D 

where D = 1 is the diameter of the hard spheres, x = rp^^^, and Xc is a selected cutoff 
distance [ISj. The rescaled radial coordinate x is set such that packing fractions of varying 
densities will be comparable, in that the total number of sphere centers over the integration 
range in packings of varying densities will be the same. The quantity Xc must be chosen 
with care: in a short-ranged disordered system, if Xc is chosen too large, the metric will not 
discriminate well between states as h{x) with x ^ oo. If Xc is chosen too small, T 
will not take contributions from correlations at greater r into account. For the purpose of 
measuring order in all of the relatively short-range pair correlation functions described in 
this paper, = 4 (or for packing fraction near 0.65, about r = 3.75) is chosen to ensure that 
all fluctuations of the total correlation functions about zero greater than a certain minimal 
amplitude (about 0.1) are taken into account. 



Using the translational order metric (19), we show here that the degree of order in the 
configurations presented is greatest for the maximum packing fraction achieved. Table |T| 
referenced by figure number and g2{r) functional form, provies the average kissing number Z 
and terminal packing fraction 0=,, versus order T for the pair correlation functions previously 
discussed. The table conveys that large differences in order may exist for similar densities, as 
is the case with the pair correlation functions represented by the top plots of Figs. [T] and |2} 
Additionally, a higher average kissing number does not necessitate higher order: the degree 
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of order in the maximum average kissing number case of the five parameter form (Fig. [3]) is 
substantially lower than in the the maximum packing fraction case (Fig. [T]). Here, greater 
correlation at a distance is sacrificed to put more particles in contact. Finally, it is of note 
that for the maximum packing fraction cases of the five- and 11-parameter forms, represented 
by the top plot of Fig. [TJ and Fig. [4], respectively, order is greatest. These results are by 
no means conclusive, but do suggest that the maximum packing fraction and maximum 
order representations of a given pair correlation function functional form are similar or even 
perhaps the same. 



TABLE I: T order metric values 
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0.631 


9.532 


0.36 




3 






13 




0.64050 


5.1593 


0.34 




2 






13 




0.64268 


5.0633 


0.43 




1 






13 




0.670 


6.672 


0.40 




6 






18 




0.6850 


10.4281 


0.46 




4 






18 





B. Moving toward crystalline order 

In the optimization study above only two sinusoids decaying as were added to the 
initial five parameter form. We now show that an infinite (or finite) number of sinusoids 
decaying as r~^, along with a Heaviside step function centered at r = 1, can represent any 
bounded (nowhere infinite) piece wise- different iable pair correlation function that decays to 
unity sufficiently fast. This follows directly from the existence of Fourier transforms, as 
will be explained shortly in greater detail, assuming that we explicitly state the bounding 
and decay conditions on the pair correlation function to mean that the Fourier transform of 
f^fi'f') = i^'^{92{i^) — 'd{r — 1)) exists. 

From a Fourier integral theorem, a radial function [37] /(r) can be represented on an 
infinite interval not containing zero, for our purposes, [l,oo), by a continuum of sinusoids 
decaying as r~^, so long as the one- dimensional Fourier transform of r^/(r) over the interval 
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exists in the sense of generalized functions [38], as follows: 



/WH-^""^'*'"*"*' (20) 

0, r < 1 



with 



1 /"°° 

F{k) = - x^f(x) cos {kx)dx. (21) 
Ji 

Under these conditions, a smooth disordered hard-sphere g2 that decays to unity can be 
written simply as 

g2{r) = f{r) + e{r-l). (22) 



For example, to represent giv{^) (Eq- (16)) in this form with C = 0, using the parameters 



from Eq. (16), we set F{k) = A6{k - B) 



Discretizing the k in expression (21) with increment Ak between two successive discrete 



values, /(r) may be approximated with F{k) a series of real constants. In the limit as 



Ak ^ or as A; becomes continuous, expression (20) holds and any (72 (r) meeting the 
aforementioned conditions may be represented by Eq. (22). Thus, to demonstrate that 
disordered packings as defined can attain densities up to tt/v^TS, it only remains to show 
that there is a realizable packing with packing fraction near vr / vTS that can be represented 



by a pair correlation function in the form of Eq. (22) for which the Fourier transform of 
r^/(r) exists, which is the focus of the remainder of this section. 

The pair correlation function of an infinite close-packed FCC crystal consisting of con- 
gruent spheres of diameter D = 1 can be written as 

1 °° f^fcc 

9tir) = —j:^S{r-q^), (23) 

where the are the distances from the origin to each successive shell in an FCC packing 
(which in these scaled coordinates are qi = y/i), the sum runs only over those i for which 
shells are present (for example, i = 14, 30, 46, etc. are skipped), the Cf'^'^ are the coefficients 
of the (even-exponent) terms in the FCC theta series pj (which represent the number of 
spheres present on a shell), and 6 represents a Dirac delta function. Any same-density 
stacking variant of an FCC configuration, i.e., any of the Barlow packings [321 SD], may be 
represented by a similar series of delta functions. For a spherical crystal of finite size, the 
remain the same but terminate for i > 21 where Vl is the distance from the central sphere 
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center to the center of the outermost sphere of the crystal, and the Cj must be reduced to 
reflect a lower average number of spheres present at given distances for all but the central 
sphere. 

An FCC configuration is not, however, disordered. To create a disordered packing from 
a series of finite-size FCC crystals, three further steps must be taken. First, the delta 
functions are replaced with smooth, sharply peaked Gaussian-like curves that decay to zero 
at distance \w/2\ from each q^. We define the integral under each of these curves from 
Qi — w/2 to qi-\-w/2 to be equal to the reduced Cj of a FCC crystal of specified finite size. 
Physically, this represents packing spheres of diameter 1 in a FCC arrangement, shrinking 
the radius of these spheres to = (1/2) — (1/4)?^', and then moving each sphere a distance 
of no more than (1/4) w from its initial position such that the pair correlation function of 
the packing near each gj has the properties just described. It is of note that the exact form 
of the pair correlation function is not of consequence as long as the proscribed properties 
are maintained. 

Second, the are again altered such that the finite-size "shrunk" FCC configuration 
of spheres becomes cubical in shape, with a cube of side length /"^''^ circumscribing the 
arrangement. Finally, an infinite number of these cubes, with principal axes of the FCC- 
packed spheres within each cube arranged at uncorrelated angles, are packed together tightly 
such that each pair of parallel faces of every cube are parallel to a pair of faces on every 
other cube, but with each face touching four other faces (i.e., four other cubes). Cubes may 
be tightly packed in this fashion such that they fill space, and such that the resulting pair 
correlation function of the spheres comprising these cubes is disordered due to the packing 
of the cubes and the random arrangement of the principal axes of the FCC-packed spheres 
within each cube. 

This pair correlation function would be in practice difficult to write explicitly due to 
boundary effects, but this is inconsequential to the description, as all that is necessary is to 
illustrate that such packings exist. Additionally, it is noteworthy that initially packing the 
spheres in each cube in arrangements equivalent to any of the Barlow packings would yield 
an equivalent result. 

Boundary effects become negligible as the number of cubes approaches infinity, and for 
increasing P"^^^ and decreasing w, any packing fraction tt/6 < cf) < n/y/l8, with n/6 — 4>^'^ 
the packing fraction of a simple cubic arrangement, can be obtained. The minimal packing 
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fraction 0**^ is obtained when only one sphere is present in each cube. The distance w can be 
made small enough for any finite-sized system of cubes such that the system exhibits physical 
stability when force is applied, though the pair correlation function of such a system would 
change (perhaps to include delta functions) as soon as stress and strain were present. 

In this way, as long as w > and P^'"^ is finite, a disordered packing may exhibit any 
packing fraction up to vr/ vTS, if one accepts that the pair correlation function of the system 



described above, represented in the form of Eq. (22), decays to unity at least as fast as r~^~*^ 
and in a form such that the Fourier transform of r'^/(r) exists. Due to this result and the 
others demonstrated in this paper, we state that the sequential addition and optimization 



of more than two sinusoidal terms of the form of expression ( 16 ) will allow packing fractions 
greater than 0.6850, and as the number of terms grows, the optimal packing fraction for a 
realizable g2 is conjectured to approach n/^/TS. 

V. CONCLUSIONS AND DISCUSSION 

Using the (yf2-invariant method with g2's satisfying the three necessary, but generally 
not sufficient, conditions for realizability, we demonstrated without implicit reliance upon 
any packing methodology that packing fractions well above 0.64 are obtainable for pair 
correlation functions incorporating the salient features of disordered packings. A packing 
fraction of 0.6850 was obtained employing a test family of (72's mimicking features observed 
in MRJ packings, including most notably core exclusion, contact pairs, and a sinusoidal 
decay to unity as r~^. Consistent with a previous study [10], we found that to achieve 
higher packing fractions, the degree of order must increase. 

Additionally we showed, employing a qualitative example and a revised version of a 
Fourier integral theorem, the surprising result that a disordered packing as defined may 
reach packing fractions approaching vr / vTS, the maximum possible for a three-dimensional 
hard-sphere packing. These results support the hypothesis that continued addition and 



subsequent optimization of sinusoids decaying as (of the form of Eq. (16)) will find 
realizable (72 with higher terminal packing fractions up to (but not reaching) vr/v^IS. 

The conclusion that the addition of sinusoids decaying as allows higher terminal 
packing fractions (where the addition of other features does not necessarily) is very relevant 
to the study of high-density physical and jammed disordered systems, since it demonstrates 
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that this feature contributes significantly in allowing the systems to reach aforementioned 
higher densities. It is noteworthy that the decay to unity present in the g2^s of MRJ 
packings also characterizes (72 's of high-density Bose systems [H], ground states of fermionic 
systems [42j , and models of the density distribution of the early Universe [l3l EH HH], Il6] , 
though these latter systems are not sphere packings. 

In future work, we will seek to extend these results to configurations of multi-component 
spheres and non-spherical objects. Additionally, we will investigate the structural origins of 
the presence of the decay in MRJ states, and attempt to determine if in higher dimensions 
additions of terms consisting of sinusoids decaying as the smallest (inverse) integer power 
for which a packing remains disordered, r^'^^^, to a g2{r) representing the salient features 
of a (i-dimensional hard sphere system, will allow for packing fractions up to the known 
maximum in that dimension. 

Moreover, we will examine the form of Zmax{R), the maximum number of sphere centers 
that can be placed within a (larger) sphere of radius R centered on an arbitrary sphere 
center. A natural question to ask is whether the realizability condition Z{R) < Z^axiR) 
further constrains g2ir) beyond the three conditions ([2]), (|3]), and (j?]). The answer is in the 
affirmative. We have already noted that the Cohn-Elkies linear programming upper bound 
formulation |22j is the dual of the Torquato-Stillinger lower bound procedure, i.e., the g2- 
invariant process |2T]. Cohn and Kumar [57] recently proved that there is no duality gap 
between the primal and dual LP programs, which means both the upper and lower bounds 
coincide when the best test functions are employed. Cohn and Elkies were able to find the 
test functions that yield the best upper bound on the maximal packing fraction in three 
dimensions: a packing fraction of about 0.778, which is well above the true maximal value. 
This means that there exists a test pair correlation function for the lower bound formulation 
that will deliver the same maximal packing fraction of 0.778, which clearly is not realizable. 
It was shown elsewhere that applying the additional condition that Z^ax must be equal 
to 12 up to some small positive a beyond contact together with the best test function in 
the upper bound brought down the maximal packing fraction appreciably from 0.778 |18]. 
This means that adding the Zmax condition to the corresponding best test pair correlation 
function will also improve the packing fraction estimate. Therefore, the Z^ax condition 
introduces additional information beyond that contained in the two standard nonnegativity 
conditions. The precise form of this additional condition has yet to be fully elucidated and 
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will be explored in future work. 
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APPENDIX A: ANALYTICAL STRUCTURE FACTOR COMPONENTS 

The analytical structure factor components, as calculated from relation (|5]), for gj, gn, 
giii, 9 IV and gv, with r = 1 the diameter of the spheres. 



Gi{k) 



An 



{k cos k — sin k) 



(Al) 



Gii(A;) 



Z 

pk 



sin k 



(A2) 




(A3) 




+ sin C[2{B + kfCi[B + k) + {B-kf{2 log {B - k) 
- log {B - kf - 2Ci{B - k))] + Ak cos k cos {B + C) 




(A4) 
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8A^/B^T „ , r. ^ /3 3 7 
15 

.5 3 9 l.o„o^^ r„ „ /I 1 5 1 



(^kB cos/c[5iF2( 



4'2'4' 



k^B^] 



+ ^^^^(4'2'4'-4^'^')]'^^"' 



(A5) 



In the expression for Giv{k), Si{k) represents the standard sine integral, Si(fc) = 
J^dx sinx/x, and Ci{k) represents the standard cosine integral, Ci{k) = — dxcosx/x. 
In the expression for Gv{k), pFg({aj}, {bj}, k) represents the standard hypergeometric pPq 
function, 

(ai)n ■ ■ ■ {ap)nk'^ 



,F,({a.},{6,},fc)=^ 



n=0 



{hi)n ■ ■ ■ {hq)nn\ 



with (oj) = 7(0 + i)/'-f{a) and 7(x) the standard gamma function. 



APPENDIX B: POWER-LAW DIVERGENCE IN NEAR CONTACT DISTRI- 
BUTION 

Evidence indicates that the addition of one of the salient features observed in MRJ-like 
states for three-dimensional hard-sphere packings, a fractional power-law divergence near 
r = 1 due to near contacts |l9l [50], does not increase packing fraction past that obtained 
from the five-parameter form. For example, the addition of the form 

gy{r) = ^-^l^e(r - 1)9(^4 - r), (Bl) 

with i?4 = 1.15 a cutoff parameter to the square-root power decay, to the gi{r), gii{r), 
giii{f) discussed earlier led to no additional increase in packing fraction under maximum 
density parameter optimization in 10,000 independent runs. Figure [7] shows a graph of the 
maximum packing fraction obtained versus the value of coefficient A4, where values plotted 
represent a random sample of 100 from the 10,000 runs conducted. One can see clearly that 
as A4 increases in value, indicative of the increased prominence of the inverse square-root 
divergence near r = 1, maximum packing fraction obtained decreases steadily. 

These results imply that the inverse half power decay near r = 1 is a feature intrinsic to 
MRJ configurations and that this feature must be diminished to increase packing fraction. 
Physically this is intuitive, as the presence of a half power decay near r = 1 indicates that 
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FIG. 7: Maximum obtained packing fraction (p versus value of parameter A^, represented by a 
random selection of 100 from 10,000 maximum density parameter optimization runs where a term 



g'y(r), given by relation (Bl), was added to the initial five-parameter pair correlation function 
functional form [13] 



there are many spheres smoothly distributed just outside of contact, i.e., that locally on 
average there is room to compress the system further. 
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